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ABSTRACT 


NOZZLE FLOW WITH VIBRATIONAL NONEQUILIBRIUM 

John Gary Landry 
Old Dominion University, 1995 
Director: Dr. John H. Heinbockel 

X, 

Flow of nitrogen gas through a converging- diverging nozzle is simulated. The flow 
is modeled using the Navier-Stokes equations that have been modified for vibra- 
tional nonequilibrium. The energy equation is replaced by two equations. One 
equation accounts for energy effects due to the translational and rotational degrees 
of freedom, and the other accounts for the affects due to the vibrational degree of 
freedom. The energy equations axe coupled by a relaxation time which measures the 
time required for the vibrational energy component to equilibrate with the trans- 
lational and rotational energy components. An improved relaxation time is used 
in this thesis. The equations are solved numerically using the Steger- Warming flux 
vector splitting method and the Implicit MacCormack method. The results show 
that uniform flow is produced outside of the boundary layer. Nonequilibrium exists 
in both the converging and diverging nozzle sections. The boundary layer region 
is characterized by a marked increase in translational-rotational temperature. The 
vibrational temperature remains frozen downstream of the nozzle, except in the 
boundary layer. 
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CHAPTER 1 


INTRODUCTION 


1.1 Background and Motivation 

Wind tunnels are perhaps the most important tools in aerothermodynamic 
research. Their usefulness depends on their ability to produce uniform flow quality. 
Wind tunnel flow is usually represented by three-dimensional flow through a high 
expansion converging-diverging nozzle. The construction of wind tunnels is costly. 
For design purposes, there is a need for computer codes which can accurately predict 
flow behavior. The prime motivation in this study will be to determine the flow 
quality of pure nitrogen gas at low pressure through a converging diverging nozzle. 

A useful property of these nozzles is the symmetry which exists in the azimuthal 
direction and allows the numerical problem to be reduced to two- dimensions. The 
fluid dynamics can be modeled using the Navier-Stokes equations. These equations 
can then be solved discretely over a finite grid. Still, with a wide variety of possi- 
ble nozzle shapes, numerical methods, boundary conditions, and representations of 
thermal properties, this problem could be modeled in many ways. Unfortunately, 
these differences may also produce much different results. Hence there is a need 
for a sound physical modeling and good n um erical technique in solving problems 
involving these flows. 


1 
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Most physical flow models axe based on the Navier-Stokes equations for fluid 
flow. The conservation of energy equation of the Navier-Stokes equations does not 
fully describe energy transfer for diatomic molecules, like nitrogen. The Navier- 
Stokes equations are used to model gases in thermal equilibrium. That is, energy 
fluctuations in the system redistribute equally am ong the degrees of freedom in the 
molecules. In an expanding nozzle, the temperature decreases rapidly as the gas 
travels past the throat. This drop is reflected in the drop in the energy associated 
with the translational and rotational degrees of freedom. The energy associated 
with the vibrational degree of freedom remains at a higher temperature that is 
closer to the stagnation temperature. Thus the system is in nonequilibrium. The 
vibrational component of the energy requires more time to lose its energy and fall to 
a temperature comparable to that associated with the translational and rotational 
energy. This finite time is the relaxation time. The ‘vibrational relaxation’ may be 
significant and therefore should be included in the physical development and the 
mathematical modeling. This is achieved by examining the energy based upon the 
translational and rotational degrees of freedom and the vibrational degree of freedom 
separately. The translational-rotational energy is assumed to be a function of a 
translational-rotational temperature. Likewise, the vibrational energy is assumed 
to be a function of a vibrational temperature. The Navier-Stokes equations are then 
modified to include an equation for the conservation of translational and rotational 
energy and an equation for the conservation of vibrational energy. 
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For equilibrium flow, the relaxation time is essentially zero and both temper- 
atures are the same. Conversely, if the relaxation time is infini te, relaxation does 
not occur and the vibrational temperature remains ‘frozen’ at some vibrational 
temperature. For a wind tunnel the relaxation time is fini te and non-zero, so the 
actual behavior lies somewhere between these two extremes [1]. The question then 
becomes, ’’What is the behavior of the vibrational temperature, and how does it 
affect the flow”? 

Many investigators have included vibrational nonequilibrium in their models 
of nozzle flow, as in references [2], [3], and [4]. Both Gnoffo [2] and Gogken [4] 
modeled ?V 2 — 0 2 gas mixtures in thermo-chemical nonequilibrium, in which they 
noted thermal equilibrium before the throat and that the vibrational temperature 
was ‘frozen’ downstream of the throat. For a test case with a very low stagna- 
tion pressure, Gogken noted vibrational freezing at the nozzle throat [4]. Also, 
the Mach number profiles displayed by Gnoffo showed reasonable unifor mi ty [2]. 
Canupp et al. noted freezing at the throat for low pressure, pure nitrogen flow 
[3]. Their comparison between equilibrium and nonequilibrium models showed that 
Mach numbers increased slightly for the nonequilibrium modeling. The flow quality 
they achieved, however, was poor. This was assumed to be due to boundary layer 
behavior. Waves were created which reflected off the nozzle wall and affected the 
flow outside of the normal ‘boundary- layer’ region [3]. Thus, the actual effect of the 
vibrational nonequilibrium on the flow quality is not clear. In order to avoid this 
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problem, flow conditions should be chosen which minimize the size of the boundary 
layer. 

These variations in results can be expected since they represent completely 
different models. Thus, their use in prediction is very limited. The only way to 
insure that the results will represent the behavior of a specific nozzle flow is to base 
one’s model on that specific flow. To ensure meaningful results, a rigorous physical 
development is also necessary. Testing low pressure flow would best display the 
effects of vibrational nonequilibrium. The Reynolds number should also be low so 
that boundary layer effects axe minim al. 

Proper modeling of vibrational nonequilibrium requires an accurate estima- 
tion of the relaxation time. The relaxation times most often used are based on 
the Millikan- White model [5], which has been shown to be accurate in the tem- 
perature range of 2000 <T< 5000° K. Unfortunately, in high expansion nozzles, 
the temperatures downstream of the throat fall well below this temperature range. 
Therefore a model is needed which describes the vibration phenomenon for low tem- 
peratures. To this end, the model of Meador et al. [6] is used, which accounts for a 
greater number of collision possibilities and a rigorous handling of second order col- 
lision effects. This improvement in relaxation time should make the nonequilibrium 
modeling more meaningful. 

Another problem in the literature is the use of the Mach number and sound 
speed in the description of the nonequilibrium flow. Mach number is calculated 
from the sound speed which is a function of the sound frequency, in nonequilibrium 
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gases. Since the sound speed is dispersive, this sound frequency in unknown. Some 
investigators use a frozen sound speed [7], which ignores the contribution of vibra- 
tional energy to the internal energy. This is only an approximation and therefore 
may undermine the integrity of the model. If the sound speed is only used to cal- 
culate Mach number contour plots to represent the results, this error may be s mall . 
However, if the sound speed is used in the numerical method [3] [7], the error may be 
significant. Also, since the behavior of the radial velocity component is dwarfed by 
the behavior of the axial velocity, contour plots of each velocity component would 
be more informative than a single Mach number plot. 

A wide variety of numerical methods are in use in solving flow problems [8]. 
The choice of a method depends on: 1) the ease of numerically modeling the physical 
problem, 2) the stability of the method, and 3) the rate of convergence to a solution. 
Most of the methods currently in use differ only in the last point. Thus, any of 
these methods would be suitable unless one required a specific amount of efficiency. 
Then one would be required to choose a method with a sufficiently large Courant- 
Friedrichs-Lewy (CFL) number. Still, no numerical method is suitable unless it 
properly models the physical equations and boundary conditions. 

1.2 Current Objectives 

The current objective is to properly simulate nozzle flow using the Navier- 
Stokes equations with the modifications for vibrational nonequilibrium. The test 
cases simulate flow with a stagnation temperature of 3000 ° K and a stagnation 
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pressure of 50 psi. The improved relaxation time is used to better model the 
nonequilibrium mechanism. 

The two-dimensional equations axe solved implicitly using both the implicit 
MacCormack method [8] and the Steger- Wanning flux vector splitting method [9]. 
Boundary conditions axe treated explicitly. The equations are solved discretely over 
an unequally spaced grid. 

The contours of the primitive variables of each case are plotted and compared. 
These results were also compared to the results obtained under the assumption of 
equilibrium flow. Special attention is given to the quality of the flow and how it is 
affected by vibrational nonequilibrium. 



CHAPTER 2 


GOVERNING EQUATIONS 


2.1 Governing Equations in General Form 

A set of governing equations are needed to relate the set of primitive flow 
variables, {/>, V,T, P}, representing density, the velocity vector, temperature, and 
pressure, respectively. For continuum gas flow, the Navier-Stokes equations of mo- 
tion are the basic mathematical model. A list of symbols appears in Appendix 
A. 

Conservation of mass: 

For a compressible gas the conservation of mass equation has the vector form: 

^+ P v-r = o, (2.i) 


where -j^ is the material derivative such that, 


D 

Dt 


| + v. v . 


Conservation of momentum: 

Each component of velocity, Vi, has a corresponding equation, 


( 2 . 2 ) 



— a 'i>j > 


( 2 . 3 ) 


7 
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where cr,y is the stress tensor, satisfying 


<7jj — PS{j T{j. 


(2.4) 


The shear stress, T{j y is a function of velocity: 


T ij “ liYiyi + Y) f,«) + ^Vk,k$iji 


(2.5) 


where rj and A are the coefficients of viscosity and 6{j is the Kronecker delta. The 
first coefficient of viscosity is calculated using Sutherland’s formula. For N 2 , 


2.16 x 10~ 8 g T? 
T + 184 


( 2 . 6 ) 


where g is the gravitational constant and T is in °R [10], and the second coefficient 
from the Stokes’ hypothesis, 

*=-§'»• (2.7) 

Conservation of Energy, one temperature: 

The change in the total energy of a fluid in motion equals the sum of the 
changes in the internal energy, kinetic energy, work done, and heat from conduction 
and radiation. For the flow investigated, radiation effects are assumed to be small 
and are neglected. The internal energy, e(T), is a function of temperature and 


de(T) 

dT 


( 2 . 8 ) 


where C v is the specific heat at constant volume. 
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The first law of thermodynamics is used to obtain the equation: 

= -PV • V + V-q-<f>, (2.9) 

where — PV • V represents the work done due to the velocity change. The second 
term, V • q, is the heat loss due to convection with 

q = —KVT (2.10) 

from Fourier’s law. The coefficient K represents the coefficient of thermal conduc- 
tivity of the gas and is assumed to be a function of temperature. The final term of 
equation (2.9) is the viscous dissipation function which has the tensor form 

* = V(T j; -V,) - V • V(r 0 ). (2.11) 

Conservation of Energy, two temperatures: 

When vibrational relaxation occurs, the energy distribution is changed as the 
system returns to a state of equilibrium. In order to model this behavior, the energy 
generated from the molecules’ vibrational degrees of freedom is treated separately 
from the energy generated from the translational and rotational degrees of freedom. 
Then, the total internal energy is written as 

e = e rt + e„, (2.12) 

which represents the sum of the translational-rotational internal energy and the 
vibrational internal energy. The molecules are assumed to be linear harmonic os- 
cillators and are excited to energy states, such that 
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where h is Planck’s constant and v is the collision frequency. Each energy state is 
populated by a density of .molecules, N{. Due to the speed of the resonnant V-V 
collisions, the populations are assumed to be Boltzm ann distributions, 


Ni = 


Ne idr 


(2.14) 


where N is the total number of excited molecules. It is also assumed, that for each 
collision, the molecules can be excited only to the next energy state. The total 
vibrational energy per unit volume can be written as 


pe„ = tiNi. (2.15) 

t 

The total rate of change in vibrational energy equals the rate of heat transfer due 
to vibration or 

§iJ v J2 Nieidv+ J s Y^ Nie M dg = l E ^r e « dv ( 2 - 16 ) 

where dv is an element of volume and ds is an element of surface area. By imple- 
menting the Gauss divergence theorem, the differential equation is obtained: 

+ E v • = E it- (2.i7) 

* t * 

where K is the average velocity of the molecules in energy state The diffusion 
velocity of a species, Ui, is defined to relate the individual velocities and the mean 
velocity, such that Ui = Vi — V. Substituting this into equation (2.17) yields 
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Using equation (2.15), the vibrational energy equation becomes 

+ V • („e„V) = -V • (£ NieiUi) + £ ^ ej . (2.19) 

* » 

The vibrational heat flux, q v , is defined by 


q v = J2 N i*iUi, 


( 2 . 20 ) 


and is substituted into equation (2.19) to give the result 


d(pe v ) 

dt 




By applying the continuity equation (2.1), equation (2.21) simplifies to 


( 2 . 21 ) 



I 


(2.22) 


The sum, ]T^ t - is the rate of change of the vibrational energy, — , which 

can be defined as 


V' &Ni ,e v (T) — e v (T v )s 

^ ei "9T = ^ 7 )- 


(2.23) 


t 

where T v is the vibrational temperature associated with e v , and T is the equilibrium 


temperature associated with e rt [6]. When T v = T, the system is considered to be 
in thermal equilibrium. The relaxation time, r, represents the time required for 
excited molecules to reach equilibrium. As with the internal energies, the specific 
heat can be written as a sum 


c v — c vr t + a 


VVI 


(2.24) 
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where C vr t and C vv are the specific heats at constant volume for the translational- 
rotational and vibrational degrees of freedom, respectively. For simplicity, define 
the quantity X such that 

/-r v — e v(T v ) 

O VV A — . 

r 

Then the final form of the vibrational energy equation can be written as 

p ^Dt^ = pCvvX ~ V ’ (2-26) 

To obtain an equation for the translational-rotational energy using the vibrational 
energy equation, equation(2.26), the elements of the one-temperature energy equa- 
tion, equation (2.9), are split into translational-rotational and vibrational compo- 
nents: 

^ £>(e r ^+e p ) + pv ^ + y ( ^+^)_ $ = ° (2.27) 

Subtracting equation (2.26) from equation (2.27) yields the translational- rotational 
energy equation: 

D(e rt ) - 

P-Jjf- + PV • V + V • q7t + pC vv X - $ = 0. (2.28) 

Each energy is related to its respective temperature by a separate specific 
heat. From the equipartition of energy principle, each degree of freedom contributes 
\ RT to the total energy [11]. Each diatomic molecule has three translational, two 
rotational, and one vibrational degree of freedom. Even for low temperatures, the 
translational and rotational modes axe fully activated, so a good approximation to 
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the specific heat at constant volume, C vr t, is § R, so that the translational-rotational 


energy has the form 


e rt — 2^' 


(2.29) 


where R is the ideal gas constant. Since the vibrational mode is not fully activated, 
the equipartition principle does not provide a good value for the vibrational specific 
heat, C vv . Substituting equations (2.13) and (2.14) into the basic definition of 
vibrational energy, equation (2.15), yields 


pe v = 


Nhv ie vtz~ 


(*+$)»» 


(2.30) 


This simplifies to 


pc v — 


NhvY.iic" 1 


ihv 


+ -Nhv. 


Ej e-m 2 

This equation can be simplified by utilizing the sums 


(2.31) 


E - 

ie — 


_ h v 

e kTv 


h v 


(1 _ e -P7V) 


(2.32) 


and 




(1 -e*rr) 


(2.33) 


Furthermore, a characteristic vibrational temperature, <f>, is defined such that 


<f> — y. Thus, equation (2.31) becomes 


N<f>ke~ & 1 ArJ , N<f>k 1 f , 

pe v = — + ~Nk<f> = — j- h -Nk<f>. 

2 - — - 


(1 - e T « ) 


e T v 1 


(2.34) 


Since Nk = pR, 


c v — 


** + In*- 


e& -1 


(2.35) 
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The specific heat at constant volume for the vibrational mode becomes 


dT v 


c„„ = ^ = R(^f 


J$L 

e T v 


Tv (e^ - — l) 2 


(2.36) 


From equation (2.25), equation (2.35) and equation (2.36) the quantity X is found 
to be 


X = 


4>r(l- e ^y 



(2.37) 


Note that when the system is in equilibrium, T„ = T, and consequentially X equals 
zero. 


The heat fluxes are functions of temperature having the form 


9rt — Krt^T, q v — K v VT v , (2.38) 

where Krt and K v are the translational-rotational and vibrational thermal conduc- 
tivities, respectively. The Prandtl number, Pr, is defined as Pr = — so 
K rt = | iv- For -^2 at 0 °C,Pr = 0.71, so K rt = 4.93rjR [12]. The vibrational heat 
flux is determined from equation (2.20). For a pure gas IV, 17, = -NDV(jj-), where 
D is the diffusion coefficient, or in this case, the self- diffusion coefficient. Then 

ql = -ND £ e.V(^) = -NDV(± £ 7V,e.) = -pDV ( i £ N iti ). (2.39) 

* t P i 

Substituting equation (2.15) into equation (2.39), 


q v = -pDVe„ = -pDC vv VT. 


(2.40) 
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The Schmidt number, Sc, is defined as Sc = and for JV 2 at 0 °C, Sc equals 0.74, 
reference [12]. Thus, from equation (2.38), 


K v = pDC v v = —C vv = 1.35t ]C vv . 

iC 


(2.41) 


The vibrational thermal conductivity is a function of both T and T v since viscosity 
is a function of T and specific heat is a function of T v . 

For N 2 , the relaxation time is taken from Meador, Williams and Miner [6]: 


where 


iW = W e 


<f> = 3395K, 6 = 3.2324 x 10 7 K,( = 95.9 K 


(2.42) 


(2.43) 


represent characteristic temperatures, and 


7(T) = /^( 1 + (1 + T )i ) 


(2.44) 


with 


±=£-)Y. 


(2.45) 


v 27^ 2 v 2 Tx J ) 

Figure 2.1 shows the relaxation time r derived from equation (2.42) as it compares 
to the relaxation time developed by Millikan- White [5]. 

Equation of State: 

The final equation is simply the ideal gas law: 


P = pRT. 


(2.46) 
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Figure 2.1. Relaxation Times of Meador et al. and Millikan- White 
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This equation does not always hold for gases at high pressures, but these departures 
from non-ideality are assumed to be negligible. 


2.2 Governing Equations in Computational Form 

The nozzle studied is defined in 3-dimensional cylindrical coordinates, (r,0,z) 
[figure 2.2]. However, there exists a symmetry in the 9 direction, which means that 
the variables do not change in the 9 direction. Thus, the ^-component of velocity, 
V$, is assumed to be zero, and the momentum equation in the 9 direction can be 
neglected. Furthermore, the derivatives with respect to 9 in the remaining equations 
are zero and can be neglected. This reduces the system of equations to six, and the 
problem to two-dimensions, (r,z). 

Symmetry about the centerline (z-axis) allows the flow to be determined from a 
single half cross-section, where 0 < r < f(z) and 0 < z < b. The function f(z) 
describes the shape of the nozzle and b denotes the length of the nozzle. The nozzle 
shape is defined in appendix B. 

To facilitate the numerical method, the each equation is written in the weak 
conservative form 


dUi 

dt 


l djrFj) dGj 
r dr dz 


+ -Hi = 0 . 


Equations (2.1), (2.3), (2.26), and (2.28) become 


(2.47) 


dp 1 djrpVr) djpVz) 

dt r dr dz 


= 0 


(2.48) 


djpVr) 

dt 


1 d 


+ " fr( r (P V r + P ~ T rr)) + ^(pV rV z ~ T rz ) - - (P + T M ) = 0 (2.49) 


dz 
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^§r + \s (r(pVrV - - T " )) + Tz (pV ‘ + p - T ") = 0 < 2 - 5 °) 

^(pe rt + j(V r 2 + V 2 )) + ~(r((p*r, + £(V r 2 + V?))W - + ,„„))+ 

• +^((/> e rt + f (V? + K 2 ))Vz - V r T„ - V,r„ + 9rU ) + = 0 (2.51) 

+ l^^ pevVr + qvr + §;( pevVz + qvt) - pCvvX = °- ( 2 - 52 ) 

The total translational-rotational energy, E r t, is defined by 

Ert = per, + |(V? + V 2 ), (2.53) 

and the total vibrational energy, E v , by 

E v = pe v . (2.54) 
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where 8 and L are characteristic lengths, Vo is the characteristic velocity, po is the 
characteristic density, and To is a characteristic temperature. The characteristic 
viscosity, 770 > is the viscosity calculated at To. The Reynolds number, Re, is defined 
as the combination Re = P° V ° L , The nozzle domain becomes 0 < r* < - and 

Vo “ 0 

0 < 2 * < The system of equations has the final scaled conservative form 

dU* 1 d{r*F*) dG * 1 rr „ 

~K7^r + ~ — jri 1 - ~ „ H H — 0 

at* r* or* oz* r* 



where U * , F* , G*, and H* are the column vectors 


u* = 


p m 

p*v; 

P *v: 

Kt 


(2.57) 


F* = 


.* 

' rr 


G* = 


L E* v J 

p*v; 

p*v r *v r * + £p* - 

P*v;v z * - ^^r* rz 

(E* rt + P*)V r * - V;^T* r - ^rV:i- t T* rz + q* rtr 
E* v V r * + q* vr 

p*v z * 

o*V*V* - — —r 

P v r v 2 Re T 1 

p*v z *v z * + p*~ £t; z 

(E^+ P-)V; - v;±r; z - + * 

E’V • + «• 

r z ' Hvz 

0 


(2.58) 


-* 

’ rz 


(2.59) 


H* = 


— Lip* 4. Ll J- T * 

gi-r + j 2 Re T ee 

0 

rLVC;.,X* 

^r~~ 


(2.60) 


rtyc:,r 
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The scaled equation of state is 


P' = p-RT'&). (2.61) 

Finally, the nozzle is mapped to the computational domain. By employing the 
change of variables, 

z* r * 

x ~ TpL' y== /C Lz^)JS' 

the nozzle is mapped to a square with domain 0 < x < 1 and 0 < y < 1, 

[figure 2.3]. The system of equations in equation (2.56) is converted to 



dU dF dG „ n 
dt* + dx + dy +H 


(2.63) 


by 


G =f(bx) 


6 F „ Lyf'(bx) c . 


f(bx) 


(2.64) 


H = ^TT G ' + -THTfC-f” + ff *)> V = V‘. 

f(bx) yf(bx ) 

When y approaches zero, y f( bx ) approaches infinity, hence the necessary limi t., 


lim 


y-o 


yf(xb) 


(. F* + H * ) = 


0 

IS fo d\ P V r ) 

« 2 JpRe dy* 

iS 6 d( T r z ) 

6* fRe dy 

I? 6V, d(r rz ) _ I} KrtqpS 2 Q*T , v 

6 2 fRe 9y 6 2 p 'dp' + Q^vv* 


L? K^q pP a 2 r„ 
6* P dy 2 


— qC vv X 


(2.65) 


is obtained from L’Hospital’s rule. 
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2.3 Boundary Conditions 

The steady solution of the equations is uniquely defined by the boundary con- 
ditions. In order to obtain a meaningful solution, boundary conditions must be 
applied that are not only applied correctly but are physically meaningful [13]. The 
conditions axe expressed in terms of the primitive variables {p,V r , V',T rt ,T.,P) 
and are converted to conservative variables when implemented numerically. 

At the inflow boundary, the translational-rotational temperature is held fixed 
at some initial value To. The vibrational temperature is assumed to be in equi- 
librium with the translational-rotational temperature and is also held constant at 
To- In computational coordinates, the gas is assumed to enter the nozzle parallel 
to the centerline, so the radial velocity is assumed to be zero. The axial velocity 
is held fixed at an initial velocity, Vo. Finally, since the flow is subsonic inflow 
and supersonic outflow, an analysis of the flow characteristics of Euler’s equations 
suggests that one boundary condition must be left free to change with the solution 
of the interior flow [14], [15]. To this end, the density values are extrapolated from 
interior data. In the computational domain, this represents the condition = 0. 
This is equivalent to the conditions |f = 0 and = 0 in the solution domain. 

At the fax-field boundary, the flow is supersonic. All variables are extrapolated, 


or 


dp _ dV z = dVr = dTrt _ dT v 
dx dx dx dx dx 


( 2 . 66 ) 



(2.67) 
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the pressure is also extrapolated. 

The symmetry of the nozzle is used to specify the conditions on the centerline. 
Assuming that the quantities above the center axis are mirrored by those below, 


dp = dV L= dV 1 = dr* dT v dP 

dr dr dr dr dr dr 


( 2 . 68 ) 


In the computational coordinates these conditions are simply 


dp dVr dV : dTrt &T V dP 

dy dy dy dy dy dy 


(2.69) 


Furthermore, along y = 0, 


V r =0 


(2.70), 


since the radial velocity at the centerline must flow equally in the positive and 
negative directions. 

Lastly, boundary conditions on the nozzle wall are governed by the nozzle 
shape and the viscosity. Due to the assumed roughness of the nozzle surface, no- 
slip conditions are applied to the velocities, 


V r = V z = 0. 


(2.71) 


The translational-rotational as well as the vibrational temperatures axe extrapolated 
such that 


dTrt = dT 1=0 

dy dy 

The pressure gradient normal to the wall is assumed to be zero, 


(2.72) 



(2.73) 
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In computational coordinates this condition becomes 

d -L = d -L 5 (!+(/'?)- = o. 

dn dx f(xb)f'(xby Ull 9y 


(2.74) 


The density on the boundary is then calculated using the ideal gas law evaluated 
on the boundary. 


2.4 Initial Conditions 

The choice of initial conditions is also important. Values that axe too far from 
the steady solution can make the numerical solution unstable. Also, the closer the 
initial conditions are to the steady solution, the less time is required for convergence. 
To minimize the convergence time, we approximate the steady state solution with 
the steady solution of the corresponding one-dimensional inviscid problem. The flow 
is assumed to be in thermal equilibrium and also isentropic. Since it is isentropic, 
there is no change in entropy, so ds = 0, where s is the entropy. Entropy is related 
to enthalpy by pTds = pdh—dP , where h is the enthalpy and h = e+RT. Therefore 
ds = 0 implies 

pdh — dP 
— — °- 

Using the ideal gas law and the fact that dh = C p dT , 

^C p dT-R^- = 0, (2.76) 

where C p is the specific heat at constant pressure. Integration of this equation leads 





± 

To 


e J o 


to 
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or 


P _ / T < t > \i/ l ~ e T 


± 

To 


JL = (-^LY( 

Po V <j>T Q ) v e £-l e 4^ _ i 


). 


(2.77) 


where To is the inlet temperature. The one-dimensional energy equation is written 


as 


C v dT + V z dV z = 0. 


(2.78). 


Integration of this equation yields 


Vz - Ko = 7iJ(To -T) + 2R<t>(—d T“ )» (2-79) 

e^o — 1 er — 1 

where V Zo is the initial axial velocity. The local one- dimensional Mach number, M, 
is calculated by the formula M = V z /a where a is the local sound speed which is 
defined a s a = (7 RT) 5 . Note that thermal equilibrium is assumed, so sound speed 
can be defined. Dividing equation (2.78) by a 2 yields 


M 2 = 


7i2T + 7 l To L)+ 7 T\%_ 


-)• 


(2.80) 


The continuity of momentum is expressed by AV z p = A*V*p *, where A is the area 
of the nozzle and the starred quantities represent the values at the throat. At the 
nozzle throat, M* = 1. FYom equation (2.77) and (2.80), the solution variables at 
the throat can be completely determined. From these values the solution variables 
throughout the nozzle are determined. Figures 2.4-2. 7 represent the solution to 
the one- dimensional model for the stagnation conditions specified for the test case 
considered in chapter 4. 




Figure 2.5. One-dimensional Model: Temperature vs Distance. 
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These values are then used as the centerline values for the two-dimensional 
viscous flow problem. They are extrapolated toward the nozzle wall. At the wall 
itself, the velocities are set equal to zero to satisfy the no-slip condition. The 
density and the translational-rotational temperature values are set to satisfy their 
specified boundary conditions. The vibrational temperature is considered to be in 
equilibrium with the translational-rotational temperature before the throat, so T v 
is set equal to T rt . After the throat, T v is assumed to be frozen and is set equal to 
the centerline T r * value at the throat. This frozen behavior is somewhat contrived, 
since the actual behavior of the vibrational temperature is unknown. Also, even 
if freezing is the actual behavior, the frozen temperature and the distance of the 
‘freezing’ point from the throat are unkn own. 



CHAPTER 3 


NUMERICAL METHODS 


3.1 Introduction 

The governing equations axe the system of non-linear partial differential equa- 
tions 


dU dF dG TT n 

ar + & + ^r +H=0 ’ 


(3.1) 


where t is understood to be t m , the scaled time and U, F, G, and H are the vec- 
tors previously defined by equations (2.57), (2.58) (2.59), (2.60), and (2.64). The 
steady-state solution to this system is desired when x and y axe restricted to the 
computational domain, 0 < x < 1 and 0 < y < 1. It is not possible to generate 
a closed, analytic solution, so a numerical solution is necessary. Many numerical 
methods attempt to accurately represent the differential system as a system of lin ear 
algebraic equations. This is achieved by replacing derivatives with fini te differences. 

For the solution variables, an initial set of values, U°, is assumed which repre- 
sent the system variables at time t. Then for some time increment, At, a solution 
set, Z7 1 , is derived from the system of equations evaluated at time t + At. The time 
dependant problem considers the discrete set of solutions, {U n : n = 0,1,2,...}, 
where U n is derived from the equations evaluated at time t + nAt. The temporal 


30 
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derivatives of the equations can be eliminated by using the finite difference: 


dU n+1 

dt 


U n+ 1 - U n 
At 


+ O(At), 


(3.2) 


which is first order accurate. By substitution into equation (3.1), an explicit method 
can be developed of the form 

U n+ 1 = F[U n ], (3.3) 


where F is some function which advances a given solution at time t + nAt to a new 
solution at time t + (n + 1)A t. Methods of the form 


U n+1 =F[U n ,U n+1 ] (3.4) 

are implicit and require different solution techniques depending upon the function 
F. 

To replace the spatial derivatives, the solution domain is discretized into a grid 
of points consisting of n + 1 rows and m + 1 columns. For every gird point (i j ) there 
is a corresponding solution set U"j. The points in rows 0 and n, or col umn s 0 and 
m, are boundary points, and the remaining points are interior points. The spacing 
between rows and columns varies to allow finer resolution near the throat and in 
the boundary layer. In the axial direction the grid is divided into four regions in 
which the spacing increases. The grid is divided into three regions in the radial 
direction, and the spacing is decreasing. This forms a total of twelve separately 
spaced regions. Using the coordinate transformation, equation (2.62), the solution 
grid is transformed into a computational grid, as in figure 3.1. From Taylor series 









Numerical Methods 33 


expansions of neighboring points, simple finite difference formulas are derived. 
Those using downstream data points are forward differences: 


dx 


Tjn — 77 " r)TI n 

^± hi dUi ' j 


Tjr 1 . 

u t,j+i 


hi 


dy 


k 1 


C/r*. 

— ^- + 0(Jfcl). 


(3.5) 


Formulas using upstream values are backward differences: 


du "i _ m - ^ 


dx 


h 2 


4- 0(h2 ), 


V?J-1 

dy k2 


4- 0(k2), 


(3.6) 


where kl, k2, hi, and h2 are the grid spacings defined in figure 3.2. 

For the Navier-Stokes equations, a popular time-dependent method is the ex- 
plicit MacCormack predictor/ corrector method which was first developed in 1969, 
reference [8]. Each time iteration requires two explicit time steps of the form: 


Predictor : 


U:t' = vt, - ^W" +1J - I?j) - g(G? J+1 - GJj) - AtH 


M 




hi 


(3.7) 


Corrector : 


l T^f 1 


5 (VZi + Vtf' - l-Wj-F"-!,)- -Gfj-i)- A^”.). (3.8) 


Explicit methods change variables locally so that information at a grid point fan 
not be allowed to travel past the next grid point in one time step. Thus a restric- 
tion must be placed upon the At time step value. This method provides a stable 
numerical solution when At satisfies the Courant-Friedrichs-Lewy (CFL) condition 
defined by 


At < 


.9(A t)cFL 
1 + 2 /Re A ’ 


(3.9) 





Numerical Methods 35 


where 


AtcFL 


-( 


Hii + Hi! 

Ax Ay 


+ o,a 


1 


+ 


(Ax) 2 (Ay) 

Ke a=m< „(f S&Eimi), 

V *7 


) 


-1 


(3.10) 

(3.11) 


and a is the local sound speed [16]. For a nozzle with large axial velocity and finely 
resolved grids, the necessary At value is unreasonably small. This results in long 
computer run times. To avoid the CFL restriction, an implicit solution method is 
used. 


3.2 Beam and Warming Method 

The most notable implicit method is that of Beam and War ming [17]. The 
temporal derivative is replaced by a more accurate, trapezoidal, time differencing 
formula: 

AU"+' = ;A1 3A ^" +1 + ^Ai^ + 0(A< 3 ), (3.12) 

where 

« 

AU n+1 = C7 n+1 — U n . (3.13) 


This allows the method to be second order accurate in time. In ter ms of AU, 
equation (3.1) becomes, 

a(AI/"+i ) - u*) _ d(F»+' - F n ) 0(G»+1 - G") , rrn+1 rrnN 

m “ m “ dx d~ y {H ~ H y (314) 


In order to simplify equation (3.14), observe that the Taylor series expansion 


of F with respect to time is 
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Since At is small, terms 0((At) 2 ) and smaller can be neglected. Also F is a function 
of the conservative variables, U, so the chain rule can be applied to to produce 


j?n+i pn dF dU 

+ dU dt 


(3.16) 


By letting A = and using equation (3.2), equation (3.16) becomes 


F n+ i =F n + AAU n+ \ 


(3.17) 


Likewise, 


G R+1 = G n + = G n + BAU n+1 . 

oU at 


(3.18) 


and 


H n+1 =H n + t = H n + CAU n+ \ 

ou ot 


(3.19) 


where B = and C = The quantities |^, |^, and are the Jacobian 
matrices. Substituting these expansions into equation (3.14) yields 

= —^(AAU n+1 ) - ^(BAU n+1 ) - CAU n+1 . (3.20) 

The derivatives in equation (3.12) can be replaced, using equations (3.1) and (3.20), 
resulting in 


AU n+1 = ^At(--^-(AAU n+1 ) - BAU n+1 ) - CAU n+1 ) + 


^At(- 


dF” dG n 

dy 


dx 


H n ). 


(3.21) 


By introducing the operator notation: 

( d A\AV»+' - (—\AU n+1 - d(* A ^ +1 ) 

K dx } dx ’ K dy } ~ dy 


(3.22) 



Numerical Methods 37 


equation (3.21) can be expressed as 


(/ + iA t§£ + 1 a t d 4 + iA«C)AI/»« = iAi(-^ 


dy 


dx 


dG" 

dy 


tf n ). (3.23) 


In the computational domain, there are mn of these equations which must be solved 
simultaneously. To achieve this requires the inversion of a block pentadiagonal 
matrix. To simplify the required computation, the operators of equation (3.23) are 
factored: 


( /+ 5 A< §7K' + iA‘f + iA<C)A^> = \At( 


dF n 

dx 


dG n 

dy 


— H n ). (3.24) 


This equation can then be solved stepwise as follows: 




(J+ + ^A<C)AC7 n+1 = AV' n+1 . 

2 dy 2 ’ 


(3.25) 

(3.26) 


Each step now requires the inversion of a block tridiagonal matrix. However, equa- 
tions (3.23) and (3.24) axe not equivalent. If the factors are multiplied, the ter ms 
4 ( A^) 2 §7 ^ + §(At) 2 |^-C arise in addition to those in the original operator. This 
produces an error proportional to ll#ll> which is insignificant for small At. 

Thus, a restriction on the At value is still required [8]. 
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3.3 Steger- Warming Flux Vector Splitting Method 

To this basic form, the flux vector splitting developed by Steger and War min g, 
reference [9], is incorporated which considers the one-dimensional Euler equations 


dU_ dF_ () 

dt dx 


(3.27) 


where F = [pV x , V x ! p + P, (e + P)V X ] T . The vector F is considered to be a homo- 
geneous function of degree one since it has the property 


F(cU) = cF(U), 


(3.28) 


for any scalar c. From the previous development, equation (3.27) can be written as 


dU_ dU_ 
dt ^ dx 


= 0 , 


(3.29) 


where A is the Jacobian matrix associated with F and which, for simplicity, is 
assumed to be constant. Since F is homogeneous it also has the property that 
F — AU . Since equation (3.29) is hyperbolic, the eigenvalues of the equation, A„, 
axe real. Also a similarity transformation exists of the form S~ l AS = E, where E is 
the diagonal eigenvalue matrix. Using this transformation and d efining U' = 5 _1 C7, 
equation (3.29) can be written as an eigenvalue problem 


dU' , x dU' 

A* A 

at ox 


= 0 . 


(3.30) 


If a backward difference is substituted for the spatial derivative, the solution of 
equation (3.30) is stable only if the corresponding eigenvalue is positive. Conversely, 
if a forward difference is used, the solution is stable only for negative eigenvalues. 
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Neither difference scheme, or any other, is stable for both positive and negative 
eigenvalues [9]. Thus directing the finite differences in the direction in which the 
information is propagated increases the stability of the system. 

This knowledge can then be applied to the present system of equations. The 
Navier-Stokes equations, however, do not exhibit purely hyperbolic behavior nor 
axe the functions F and G homogeneous. To this end the functions F, and G axe 
split into homogeneous functions, Fh and Gh, and non- homogeneous functions, F v 
and G v , such that 


F — Fh + F v , G = Gh + G v , 


(3.31) 


and 


F ‘ = l 


r* 

z 


r * 
' z 


p*v> 

P *v;v 2 
p*v;v; + p 
(E; t + p*)v : 


F — ^ 
Fv ~J 


> Re T rz v z Re T zz + Hr tz 


E*V* 

v z 

0 

6 2 Re T rz 

L 2 1 _* 
~WTu t zz 

* 1 


Gh = 


f(ib) 


L 

+i’„ 

J 

p*v; 


‘ p*v; - 

p*v;v; + 


p*v;v; 

p*v;v; 

Lyf(xb) 

K lb ) 

p*v;v; + p* 

(E; t + P')v; 


(eu + p*)v: 

E' v v; 


- e*v; . 


(3.32) 


(3.33) 


(3.34) 



Numerical Methods 40 


0 


G v = 


f(xb) 


-L 2 -Lt* 

W Re T rr 

iJ_r* 

Re T rz 

r Re'rr 6 2 V z R t T rz 


+ <l*rtr 


< r 


Lyfjxb) 

f(xb) 


0 

L 2 1 _♦ 

~W7U T rz 

L r * 

Re T zz 

— V*-i-r* -4-n* 
K r Re'rz v z R e T zz ' Vrt 


(3.35) 


Also define Ah = Bh = A„ = , and note that A — Ah+A v 

and B = Bh + B v . The vector H and matrix C are not altered. 

To implement the vector splitting for matrix Ah, the eigenvalues of Ah are 
calculated numerically. From these the diagonal eigenvalue matrix E is formed. 
Also, the matrix S and its inverse S -1 of the similarity transformation are calculated 
such that Ah = SES -1 . Next, the diagonal matrix E is split by E = E + +E~, where 
E + contains the positive eigenvalues and E~ contains the negative eigenvalues. 
Finally, the flux matrix Ah is split by Ah = A + + A~ where 


A + = SE + S~\ A~ = SE~S~ 1 . 


(3.36) 


Similarly, the flux matrix Bh is split using its respective eigenvalues such that 
Bh = B + + B~ . This allows equation (3.24) to be written as 
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■('+ + 9 -f + &»<'+ l^ a S- + ^ + ^ > + ^> Atr+i 


1 A ±( dF dG tt \ | n 


(3.37) 


This is equivalent to 


(' + 5 Ai (^ ' + d -& + £))*'*” = ;*(-£ - f - ff) r. (3-38) 


< 7 + s^tt + fr + + i A,c > At/ " + ' = AV " +1 - 


(3.39) 


Since A + and B + represent the homogeneous contributions of the system with 
positive eigenvalues, backward differences are used for and Likewise, 

forward differences are used for and ^JL-. The non-homogeneous matrices A v 
and B v are not affected by flux splitting. They are primarily viscous terms and 
central differences are used for and since they do not create dispersive 
instabilities [14]. For the given grid structure, the central differences are 


dA v ij 1 { h2 A ,hl h2, A hi . 

'[rr^+bi ■*" (^2 — ~ To^-vi- i,j] + 0(h ), (3.40) 


dx 


hi + h2 L hl 


h 2 


and 


dB 


Vt,J 


dy 


kl + k2 l k2 


r fc2_ ,fcl k2, n kl _ , _ /t ov . . 

[~jP)Bvi,j+i + ( j- 2 — ~ )> (3-41) 


k2 


where h = max(hl, h2) and k = max(kl, k2). Since the errors associated with these 
formulae are quadratic, they provide more accuracy than the forward or backward 
differences. For this reason, spatial derivatives that are evaluated in the Jacobian 
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matrices as well as the F, G, and H matrices are also central differenced when 
possible. By replacing the derivatives in equation (3.38) and equation (3.39) by the 
appropriate difference formulae, the equation for becomes 


AV^ 1 + 


+ 


J A ‘(j5 (A & A, tf' 1 “ AV '-" )+ Tk (A *< AV Sn. ~ WKt'* 


A£, y 1 + i a, (s (a 6 A0 « 1 - *&-i AO tf- 1 >> + ^( B r,-. At ?j"-. - »,^C)+ 

+ (| - {£)**, -]) - 


1 rk2 


kl -f- k2 


Rearranging terms yields 


( 3 . 43 ) 


h 2 


1 , 1 

2 'hL i+1 ’i + hl(hl + h2) 


A„ l+lii )AV^. + 


+ (l + -At(— A + L^- j_ * ^2 

+ ^ + 2 ^2^*-' + 

1 . ./ 1 hi 


AV’I'^j = rhs( 3.38), (3.44) 


+ 2 A< ^ h2 A ^- 1 ’’ h2(hl + /i2) 


1 1 A;2 

2^M Bij+1 + tl(H + k2) Bm ' ,+ ^ AU '‘l+ 1 + 

+(! + - |)^))a^-+ 

11 itl 

+ 2 Ai (-it2 B i-' - S J Jr ^lW..)AD!S? 1 = (3.45) 

This system of equations has the matrix form: 


[C,][AV" +1 ] = (rH. 


(3.46) 


[C 2 ][At7" +1 ] = [AV n+1 ), 


(3.47) 



Numerical Methods 43 


where [C\] and [C 2 ] are (5nm) square coefficient matrices and [AV n+1 ], [AI/ n+1 ], 
and [rhs] axe column vectors of length 5nm. The matrices [Ci] and [C 2 ] axe block 
tri diagonal and can be readily inverted. Figure 3.3 shows the non-zero entries of 
the matricies [C\\ and [C 2 ], where each element is a 5x5 matrix. 

For some values of i and j, the AV rn+1 and A U n+1 values associated with 
boundary points axe required. In cases where the boundary values remain constant, 
AV n+1 = 0 and AU n+1 = 0, and the form of equations (3.44) and (3.45) do not 
change. However, for complicated boundary conditions, the form can be signifi- 
cantly changed so that [C\\ and/or [C 2 ] are no longer tridiagonal. To avoid this 
problem, the boundary points axe treated explicitly. In the implicit step, backward 
or forward differences axe chosen at near boundary points so that the boundary 
points axe not required. This preserves the tridiagonal nature of the matricies [Cl] 
and [C2]. The boundary points axe only used in the explicit step. Since every 
time step is restricted by At, the variable changes at one time step are assumed to 
be nearly equal to the changes at the previous time step, so AV”"*" 1 « AV n and 
AU n+1 ~ A U n . Thus, boundary value changes can be calculated from known infor- 
mation and can be moved to the right hand sides of equations (3.46) and (3.47). To 
insure that the boundary conditions axe satisfied, the boundary values are updated 
after each solution of the system. 

The system can then be solved numerically and the set of AU n+1 values ob- 
tained. Using 


U n+1 = U n + AU n+1 , 


(3.48) 
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Cl,l Ci ,2 

^ 2,1 C 2,2 <^ 2,3 


^ 3,2 ^ 3,3 073,4 


07*4,3 £4,4 0*4,5 


0*5,4 05,5 0*5,6 


0*6,5 0*6,6 0*6,7 

0*7,6 0*7,7 0*7, J 


0*8,7 08,8 08,9 


L 09,8 09,9 

Figure 3.3. a) Pattern of Non-zero Elements of Matrix [Cl]. 

b) Pattern of Non-zero Elements of Matrix [02]. 
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the variable values at the next time step are found by updating the previous values. 

3.4 Implicit MacCormack Method 

Another implicit method that is easily implemented is the implicit MacCormack 
method. As in the other implicit methods, the Taylor expansion of the governing 
equations resulted in equation (3.20), Using equation (3.1), this equation can be 
written as 

ft B ft F n ftG n 

AU n+l +A*( — (AAU n + 1 ) + — (BAU n+1 ) + CAU n+1 ) = At ( — H n ). (3.80) 

ox oy ox oy 


This equation is then solved in the following sequence: Predictor step: 

(Alt ij) B ,*ci. = -Ai(^T + + (3 ' 81) 

(. I + + Ai^A + A(C)Att*” +1 = (AC/") E , pl , c ,„ (3.82) 

U*p +1 = UPj + At/* n+1 (3.83) 

Corrector step: 

(AUTjh.pH'i, = -A + V ^ );i + 0 a * )5) (3.84) 

(I + + A <^^ + A ‘C)AC'," +1 = (At 7j)) Explici , (3.85) 

or* 1 = i[tt; + (o%*'av?*'] (3.86) 


Here A x and A y represent forward differences as in equation (3.5), and V x and V y 
are backward differences similar to equation (3.6). 
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The solution of these equations requires solving block tridiagonal matricies 
which are either upper or lower tri diagonal [figure 3.4]. This, however, causes no 
problem in implementation. The boundary points are treated explicitly. Since every 
time step is restricted by At, the variable changes at one time step are assumed 
to be nearly equal to the changes at the previous time step, so AU n+1 ~ A U n . 
Thus, boundary value changes can be calculated from known infromation and can 
be moved to the right hand sides of equations (3.82) and (3.85). To insure that 
the boundary conditions are satisfied, the boundary values are updated after each 
solution of the system. 

The major differences between this method and the Steger- Warming method 
is that the Jacobian matricies are not split, so convergence is assumed to be slower. 
In both methods, each iteration requires two ‘sweeps’. The implicit MacCormack 
method sweeps in the downstream direction and then the upstream direction. The 
Steger- Warming method sweeps in the x-direction and then the y-direction. This 
difference may have some effect on the results. 

3.5 Numerical Implementation of Boundary Conditions 

To insure that the boundary conditions are satisfied, the boundary point val- 
ues are updated after every iteration. The conditions specified in chapter two are 
converted from differential to algebraic equations using finite differencing. 

Along the centerline, j=0, the conditions are of the form ^ = 0, which are 
equivalent to ^ = 0. This is implemented by replacing ^ by a simple forward 
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Figure 3.4. a) Pattern of Non-zero Elements of Leading Matrix in Predictor Step. 


b) Pattern of Non-zero Elements of Leading Matrix in Corrector Step. 
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difference 

+ = °' (3.87) 

The numerical boundary condition then becomes 

U it0 = Ui,i + O((kl) 2 ). (3.88) 

Due to symmetry, V r = 0, so U 2 ,i,o is simply set equal to zero. 

Along the far field boundary, i=m, the extrapolation conditions are of the form 
= 0. A backward difference is used to formulate the numerical condition 

U m ,j = U m - i,j + 0((h2) 2 ). (3.89) 

At the inflow boundary, i=0, all of the variables are kept at their initial values 
except the density which is extrapolated using the linear extrapolation similar to 
that used on the far field boundary, 

Ui,oj — ^i,i, j + 0((hl) 2 ). (3.90) 

On the nozzle wall boundary, j=n, the no-slip conditions are enforced and 

U 2 ,i,n = U 3 ,i <n = 0. (3.91) 

The temperatures at the wall are extrapolated linearly, 

T.,„ = + 0((H) 2 ), 


(3.92) 



Numerical Methods 49 


where T is either Trt or T v . Using equation (2.74), the wall condition, = 0, can 
be written as 


dP _ 6(1 + (f'(xb)) 2 ) dP dP 


= a- 


dx f(xb)f'(xb) dy dy ' 


Replacing with the second order backward difference 


3 Pjj ^Pj,j-i + Pj,j - 2 

2k2 


+ 0((k 2) 2 ), 


and ^ by the backward difference 


~ +0(h2), 


equation (3.93) yields the munerical condition 

P itj =a -1 +a(4P, tJ _i - P itj _ 2 )] +0(max( 


k2(h2) 2 (k2) 3 h2 

k2 + h2’ k2 + h2 


(3.93) 


(3.94) 


))• (3.95) 


Since k2 is usually smaller than h2 at the wall, the error term is at most 0((62) 2 ), 
therefore it is of the same order as the rest. The density can then be calculated 
from the wall pressure and temperature using the ideal gas equation. 


3.6 Conclusion 

These methods require a certain amount of bookkeeping to be implemented. 
They are second order accurate in time and first order in space [16]. The use of the 
implicit method removes the CFL restriction on the time step. However, restrictions 
on the time step and grid spacings are still required, to keep the computational error 
small. Flux vector splitting is added to increase the stability of the munerical 
method. The boundary points are handled explicitly in order to maintain the 
structure of the matrix equations to be solved. 


CHAPTER 4 


NUMERICAL RESULTS 

4.1 Introduction 

To discover the effects of thermal nonequilibrium on nozzle flow, a test case 
is solved numerically using the methods discussed in the previous chapter. The 
test case had a stagnation pressure of 50 psi and a stagnation temperature of 3000 
degrees K. Since there is assumed to be an initial velocity at the entrance of the 
nozzle, the term ‘stagnation condition’ is a misnomer. However, this velocity is 
comparatively small to any characteristic velocity, so these initial conditions are 
nearly the stagnation conditions. Also, some preliminary considerations must be 
investigated. 

4.2 Validity of the Navier Stokes Equations 

Rapidly expanding flows are characterized by a significant decrease in the den- 
sity of the gas beyond the throat. For flows with low stagnation pressures, the 
exit densities can be extremely low. If the gas is very rarefied, it may no longer 
be a continuous fluid. In this case, the Navier-Stokes equations would no longer 
accurately represent the flow. In that case, the rarefied gas flow would usually be 
solved using Monte-Carlo techniques. 


50 
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In order to determine the amount of continuity in the gas, the Knudsen number, 
K n , is employed. The Knudsen number is a dimensionless quantity defined as 
K n = •£, where A is a mean free path and L is a characteristic length of the flow. 
For nitrogen gas, the maximum mean free path is approximately 



where <7 2 = 14.9(10) -16 cm 2 and n is the density of the gas in molecules / cm 3 [11], 
The maximum mean free path of nitrogen versus the gas density is represented in 
figure 4.1. The choice of a characteristic length is nonspecific, so the radius of the 
nozzle is selected. Since the nozzle radius is not a constant, the local Knudsen 
number at some distance z from the nozzle entrance is defined as the mean free 
path evaluated at the centerline point at z divided by the radius of the nozzle at 
z. From the one- dimensional model, the mean free path can be estimated for the 
two-dimensional test case. The local Knudsen number is calculated, and is plotted 
for all points throughout the nozzle in figure 4.2. 

According to Bird, [18], the Navier-Stokes equations can be used for flows which 
have Knudsen numbers of 0.1 or lower. For flows with greater Knudsen numbers, 
the Navier-Stokes equations break down and do not accurately represent the flow. 
Therefore, there is a lower limit as to what inlet density can be used. From figure 
4.2 it can be seen that the test case lies within this region of usability, and the low 
stagnation pressure should not cause any significant problem. At most, the results 
obtained at the nozzle exit plane would be suspect. 



Maximum Mean Free Path [m] 




Figure 4.2. Local Knudsen N 
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4.3 Boundary Layer and Turbulence 

Another consideration is the size of the boundary layer and the existence of 
turbulent behavior in the boundary layer. Boundary layers exist for viscous flows 
but the width of this layer is a function of the flow characteristics. The width, 
can be approximated by 


X , c ^ c 


(4.2) 


pV' (Jfc)*’ 

where C is a characteristic length [19]. Using the same argum en t as for Knudsen 
number, a local boundary layer width can be found as a function of position in the 
z-direction. As before, the radius of the nozzle is used for the characteristic length. 
Density, velocity, and viscosity values are evaluated at the centerline. Figure 4.3 
represents the size of the boundary layer at all points along the nozzle. From this it 
can be seen that the width of the boundary layer is relatively small. Therefore, be- 
havior in this region should be relatively isolated and should not have a pronounced 
effect on the majority of the flow. 

Since turbulent behavior occurs in the boundary layer, turbulence modeling is 
not included in the model. The benefit of increased accuracy in the description of 
the flow, in such a small region, is not great enough to warrant the added complexity 
of the numerical model. However, it should be noted that Korte et al. found that 
not resolving the behavior at the throat affected the downstream Mach number 
(velocity) values [20]. Properly resolving this region, however, requires using a grid 
structure at the throat which is finer than usually used. Since the grid structure 



Boundary layer width [m] 
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Distance from inlet [m] 


Figure 4.3. Boundary Layer Width vs Distance from Inlet. 


Numerical Results 56 


in figure 3.1 has the greatest resolution at the throat, this may limit the growth of 
this problem. 

Since the region outside the ‘turbulent region’ is usually treated as non-turbulent , 
and only matched to the turbulent region, it can be ass um ed that the overall effect 
on the near-centerline values is minimal. In terms of wind tunnel flow, it is not as 
much the behavior of the turbulent region that is of interest, but the spatial limits 
of the section with uniform flow. 

4.4 Numerical Solutions 

The solution of the test flow was determined for three separate cases. The first 
case was solved on a 81 x 41 grid (3321 points) using the implicit MacCormack 
method. Case two used a 81 x 41 grid and the Steger- Warming flux vector splitting 
method. The third case also used vector splitting but used a different wall boundary 
condition for the vibrational temperature. These three cases were compared with 
each other and with the solution of the one temperature, equilibrium model. 

The first case, using the implicit MacCormack method, was run for 7000 it- 
erations using a scaled time step of 1.0 x 10 -3 . This translates to a total time of 
7.0 x 10 -7 s. Figures 4.4 through 4.15 are the set of contours for the set of six 
primitive variables. They also include enhancements of the throat region. Values 
stated below are centerline values. 

In figures 4.4 and 4.5, the density is seen to fall off greatly from .373 km/m 3 to 
3.26 x 10 -5 fcm/m 3 . This mirrors the pressure drop of 3.3 x 10 5 Pa (50 psi) to 
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1.08 Pa [figures 4.6 and 4.7]. The ratio of exit pressure to initial pressure is small 
enough so that a shock should not occur between the throat and the exit plane. 

The radial velocity contours, figures 4.8 and 4.9, match intuitive expectations. 
The V r values are negative before the throat (converging section) and positive af- 
ter the throat (diverging section). The values are' considerably less in magnitude 
compared to the axial velocity values, so they would not affect the vector value of 
total velocity. The absence of turbulence is also noticeable. This, however, is due 
to restrictions placed upon V r in the numerical routine. Since the AU n+1 value 
calculated in the solution routine for V r usually called for at least a 10 6 percent 
change in the V r value, instabilities would arise not only in the boundary layer but 
also outside the boundary layer. Since the V r values were relatively small, the need 
for numerical stability outweighed resolving the exact behavior of the V r contours. 
Therefore, the V r values were not allowed to change more than one hundred percent 
per iteration. This, of course, would limit the growth of turbulence. 

The axial velocity contours, figures 4.10 and 4.11, showed that the velocity 
increased from 6 m/s to 2837 m/s. The flow is uniform with almost linear contours 
throughout most of the nozzle, and then rapidly decreasing in the boundary layer 
to satisfy the no-slip condition. The size of the boundary layer is small as expected, 
and causes no apparent problems. 

The translational-rotational temperature falls sharply after the throat which is 
expected since the pressure is also dropping [figures 4.12 and 4.13]. At the far field 
boundary, the exit temperature is 112.07° K. This, along with the low pressure, 
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still represents gaseous nitrogen. Boundary behavior is very noticeable. It is the 
same size as the boundary layer seen the V z contour. The temperature is increasing 
to the wall. To understand why this occurs, note the equilibrium energy equation 


p-z - + pV • Ve + -PV - V + V- 9 + $ = 0! 
a t 


For the values at the wall, the no-slip condition makes V r = V z = 0, so the viscous 
dissipation terms and pressure drop terms are zero. Equation (4.3) then becomes 


de 


Since V • q = — V • (KVT), where K is the thermal conductivity, equation (4.4) can 


be written as 


= V • (KVT). 


If it is assumed that the equilibrium temperature, T, decreases at the wall, then 
VT is negative, so the internal energy would decrease. This would imply that 
the temperature would also decrease. However, the translational-rotational energy 
equation from the nonequilibrium model at the wall boundary, differs due to the 


coupling term, 


r\ 

P—fif = v ’ (KVT) — pC vv X. 


This can be written as 


V • (KVT) - p( 6p(T) r et,(Tp) ). 
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Since T <T V and T > 0, the difference in the last term of equation (4.7) is negative, 
so 

= V • (KVT) + p\C„X\. (4.8) 

The coupling term increases the internal energy. If this term dominated the heat 
loss term, then one should see an increase in the internal energy and, therefore, an 
increase in the translational-rotational temperature. This would cause the gradient 
of T in the r direction to be positive and V-(ATVT) would not necessarily be negative 
or at least would be less negative. Thus, numerically, one effect would create and 
encourage the other. Physically, this effect is still meaningful. As the gas slows at 
the wall, there is more time for the gas to relax. More relaxation translates into a 
greater transfer of vibrational to translational-rotational energy. This would cause 
an increase in the translational-rotational temperature. 

The most notable feature of the vibrational temperature contours, figures 4.14 
and 4.15, is the freezing of the vibrational temperature downstream of the throat. 
Also, the flow is not in equilibrium before the throat, as other authors have suggested 
reference [2]. The translational-rotational temperature drops much faster in the 
converging section than the vibrational temperature. Therefore, the assumption of 
equilibrium flow before the throat should not be used. In the boundary layer, the 
vibrational temperature decreases at the wall up to the ‘freezing point’. After this 
point, the temperature freezes and the boundary behavior disappears. The decrease 
occurs due to the slow speed, which allows the vibrational temperature to relax and 
move closer to the translational-rotational temperature. The fact that this behavior 
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stops farther downstream may mean that the translational- rotational temperature 
is increasing faster than the vibrational temperature is decreasing. 

Case two produced similar results. Using the Steger- Warming flux vector split- 
ting, the shape and the values of the contours were nearly the same. The densities, 
figure 4.16, dropped from 0.371kg /m 3 to 3.27 x 10 ~ 5 kg/m 3 . The pressure dropped 
from 3.3 x 10 5 Pa to 1.08 Pa [figure 4.17]. 

The radial and axial velocity contours, figures 4.18 and 4.19, also were nearly 
identical. The same restriction was applied to the V r changes as in case one. The 
axial velocity increased to 2836 m/s, which differed case one by only one m/s. 

The translational-rotational temperature [figure 4.20] dropped to 112.2 ° K. 
The vibrational temperature behaved similarly to case one [figure 4.21]. The vibra- 
tional temperature was in nonequilibrium before the throat and freezing occurred 
downstream of the throat as in case one. The freezing temperature was exactly 
the same as in case one. The only difference was the point of freezing, which was 
slightly further downstream. 

Overall, the behavior was similar enough to validate the results in case one. 
This case also produced uniform flow. The freezing of the vibrational tempera- 
ture in the boundary layer downstream of the throat reaffirms the idea that the 
translational-rotational temperature is increasing faster than the vibrational tem- 
perature is decreasing. 

To remove the problems encountered in cases one and two, a third case is exam- 
ined using different boundary conditions. The stagnation conditions and boundary 
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conditions were the same as the previous cases except at the wall boundary. At this 
boundary, only the condition on the vibrational temperature was changed. Since 
the gas is not moving at the wall, the energy must be in equilibrium. Therefore, 
the translational-rotational temperature and the vibrational temperature must be 
the same. To enforce this statement, the vibrational temperature is explicitly set 
equal to the translational-rotational temperature. 

The contours for density [figure 4.22] showed that the density dropped from 
.371 kg/m 3 to 3.25 x 10 ~ 5 kg/m 3 . The pressure [figure 4.23] dropped equivalently 
from 3.3 x 10 s Pa to 1.08 Pa. 

A main difference in this case was the radial velocity [figure 4.24]. which showed 
some turbulent behavior. This turbulence begins at the nozzle throat within the 
boundary layer region. Downstream the effects axe less noticeable. This is the 
behavior that was expected for the boundary layer region. 

The axial velocity [figure 4.25] contours were similar to cases one and two. The 
velocities rise steadily from 6 m/s to 2850 m/s, which was higher than the other 
cases. The contours axe uniform throughout the nozzle. 

The translational-rotational temperature [figure 4.26] drops from 3000 ° K to 
112.1 0 AT. In the boundary layer region, the temperature still increases, but much 
less severely. The vibrational temperature [figure 4.27] freezes at the same frozen 
temperature as in cases one and two, but the boundary layer behavior is different. 
The vibrational temperature drops to the value of the translational-rotational tem- 
perature on the wall boundary. Thus, the boundary behavior which stops at the 



Numerical Results 62 


throat for the other cases, now continues to the nozzles exit plane. This is the re- 
sult that was expected. Therefore, the change in boundary conditions significantly 
improved the results. 

For comparison purposes, the results of the three cases were compared to the 
results of a nitrogen flow with the same stagnation conditions but with the unrealis- 
tic assumption of forced thermal equilibrium [21]. Figure 4.28 shows the difference 
between the centerline values for the equilibrium temperature compared to the cen- 
terline values for the two temperatures found in case one. The actual measurements 
of T v are difficult in a laboratory situation, since the gas would equilibrate upon 
contact with a temperature sensor. Therefore, it would be useful to convert the 
translational-rotational temperature and the vibrational temperature into a single 
temperature that could actually be measured. Park suggested a geometric average 
= T avg , where 0 < a < 1, in order to calculate reaction rates [22]. This is 
just a rough estimate, so it is not used here. 

Other differences between case one and the equilibrium model results were 
the absence of boundary layer behavior for temperature, and an increased axial 
velocity for the equilibrium model. In this model, the exit velocity was 2900 m/s. 
This increased value may have been due to an improperly imposed conservation of 
mass flow condition in that program. That model also predicted uniform flow. 
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Figure 4.4. Case One: Logarithm of Density Contours. 
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Figure 4.5. Case One: Logarithm of Density Contours (throat region). 



Numerical Results 65 


r 


0.020 


0.010 


0.000 


- 0.010 


- 0.020 



P [Pa] 





330333 

313816 

297300 

280783 

264267 

247750 

231233 

214717 

198200 

181684 

165167 

148650 

132134 

115617 

99101 

82584 

66067 

49551 

33034 

16518 

834 

75 

17 

5 

3 

1 


0.0 0.1 0.2 0.3 

Z 


Figure 4.6. Case One: Pressure Contours. 
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Figure 4.7. Case One: Pressure Contours (throat region). 
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Figure 4.8. Case One: Radial Velocity Contours 
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Figure 4.9. Case One: Radial Velocity Contours (throat region) 
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Figure 4.10. Case One: Axial Velocity Contours. 
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Figure 4.11. Case One: Axial Velocity Contours (throat region). 
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Figure 4.12. Case One: Translational-rotational Temperature Contours. 
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Figure 4.13. Case One: Translational- rotational Temperature Contours (throat region) 
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Figure 4.14. Case One: Vibrational Temperature Contours. 
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Figure 4.15. Case One: Vibrational Temperature Contours (throat region). 
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Figure 4.16. Case Two: Logarithm of Density Contours. 
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Figure 4.17. Case Two: Pressure Contours. 
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Figure 4.18. Case Two: Ra-dial Velocity Contours. 
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Figure 4.19. Case Two: Axial Velocity Contours 
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Figure 4.20. Case Two: Translational-rotational Temperature Contours. 
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Figure 4.21. Case Two: Vibrational Temperature Contours (throat region) 
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Figure 4.22. Case Three: Logarithm of Density Contours. 
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Figure 4.23. Case Three: Pressure Contours. 
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Figure 4.24. Case Three: Radial Velocity Contours (throat region) 
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Figure 4.25. Case Three: Axial Velocity Contours 
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Figure 4.26. Case Three: Translational-rotational Temperature Contours 
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Figure 4.27. Case Three: Vibrational Temperature Contours. 
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Figure 4.28. Centerline Equilibrium and Nonequilibrium Temperatures. 


CHAPTER 5 


CONCLUSIONS 


The modeling of nozzle flow is a complicated problem. A wide variety of pos- 
sible conditions and assumptions lead to a variety of flow characteristics, which 
makes general comparisons and prediction less meaningful. However, for specific 
cases, computer simulations are cost-effective tools for nozzle and wind t unn el re- 
search. 

The governing equations for fluid flow were derived for a gas in thermal nonequi- 
librium. The internal energy is split into the energy from the translational and 
rotational degrees of freedom and the energy from the vibrational degree of free- 
dom. Each energy is assumed to be a function of a separate temperature. The 
two-dimensional viscous Navier-Stokes equations were modified to include two en- 
ergy equations coupled by means of the relaxation time. The relaxation time was 
derived from a new expression which was more realistic for the low temperatures 
that are characteristic of high-expansion nozzles. The other thermal properties were 
defined to be physically realistic and a function of temperature. 

The continuous system of nonlinear partial differential equations were replaced 
by difference equations and were solved over a discrete grid. The implicit 
MacCormack method and the Steger- Warming flux vector splitting method were 
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applied to the equations. A set of boundary conditions were used which were 
physically meaningful, and gave a unique solution. 

To test the effects of vibrational nonequilibrium, the solution of low pressure 
nitrogen gas flow through the given nozzle was tested. A stagnation pressure of 
50 psi and stagnation temperature of 3000 ° K defined the test case. Results were 
found using both methods. The results showed that unif orm flow existed for all 
cases, except in the boundary layer region. The axial velocities increased from 6 
m/s to approximately 2837 m/s. Translational-rotational temperatures fell from 
3000 ° K to 112 ° AT, except in the boundary layer which saw a marked increase in 
temperature. This was due to the addition of the coupling term to the translational- 
rotational energy equation. Properly defining the vibrational temperature boundary 
conditions reduced the negative numerical results. 

A low pressure nitrogen wind tunnel would be a viable research tool, as it should 
produce mainly uniform flow. Since turbulent behavior was restricted, the actual 
behavior, and the extent to which it disrupts the flow is still not known. Based 
upon the relatively small radial velocities and the predicted size of the boundary 
layer, it is doubtful that turbulence would have had a pronounced effect on the flow. 

The translational-rotational temperature falls throughout the nozzle, especially 
after the throat. In the boundary layer, the translational-rotational temperature 
increased significantly. The surplus was due to the vibrational relaxation, which is 
more complete in the boundary layer. Proper boundary conditions must be used 
to insure that the temperature does not increase unrealistically. Also, some of the 
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surplus energy may be due to turbulence. Proper turbulence modeling may have 
accounted for some of the excess energy and limited the growth of this temperature. 

In the diverging section of the nozzle, the gas was found to be in nonequilib- 
rium. The vibrational temperature remains higher than the translational-rotational 
temperature. This means that the assumption of thermal equilibrium before the 
throat is incorrect. 

Downstream of the throat, the vibrational temperature freezes. This ‘frozen’ 
behavior extends from the centerline to the boundary layer. In the boundary layer, 
the vibrational temperature decreases to the value of the translational-rotational 
temperature. 

A point in question is the formulation of the vibrational model. It is somewhat 
questionable as to whether the thermal properties can be split into two components 
as easily as was done here. Also, one would expect that the equation of state would 
be a function of the two temperatures. Finally, a reasonable relationship between 
the two temperatures and a physically meaningful single temperature would be 
useful for comparison purposes. These questions are areas for further research. 
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APPENDIX A 


LIST OF SYMBOLS 


a 

A 

b 

C 


Cxi i Cxjrt i Ct 


D 


D 

Dt 


d - 

« +v 


^rt •> 


f 

h 

h 

hi, h2, kl, k2 

k 

k,k„,k v 


Kn 


L,6 


Speed of Sound [ m/s ] 

Cross sectional area [m 2 ] 

Nozzle length [m] 

Characteristic length [m] 

Specific heat at constant volume [J/KgK] 
Diffusion constant 
Material or substantial derivative 
Energy per unit mass [J/Kg] 
function: radius of nozzle [m] 

Planck’s constant [ J s] 

Enthalpy [J/m 3 s] 

Irregular grid spacings 
Boltzmann’s constant [J/K] 

Thermal conductivities [IF/m K] 

Knudsen number 
Scaling lengths [m] 
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m 

M 

n 


Ni 


Qi Qrt 5 Qv 


p 


Pr 

R 

Re 
r , z 

s 

Sc 

t 


T,Trt,T v 

V 

V r , V e , V 2 

W 

X 


*,y 


Number of grid points in axial direction 
Mach number 

Number of grid points in radial direction 
Population densities [1/m 3 ] 

Heat input per unit volume [J/m 3 s] 
Pressure [N/m 2 ] 

Prandtl number 

Gas Constant [J/KgK] 

Reynolds number 

Physical coordinates 

Entropy per unit volume [J/m 3 s K ] 

Schmidt number 

Time [s] 

Temperatures [K] 

Velocity [m/s] 

Velocity components [m/s] 

Molecular weight of N 2 [Kg] 

Coupling term [K/s] 

Computational coordinates 
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V 

A 


''max 

e 

T 

hi 

A*, Ay 

Vx,Vy 


e; 


(Tij 


‘V 


Mt 


$ 



Viscosity coefficient [Kg /ms] 

Second coefficient of viscosity [Kg /ms] 
Maximum mean free path [m] 

Density [ Kg/m ; 3 ] 

Relaxation time [s] 

Boundary layer width [m] 

Forward difference operators 
Backward difference operators 
Energy states [ J] 

Stress tensor [ N/m 2 ] 

Shear stress tensor [N/m 2 ] 
Characteristic temperatures [K] 

Viscous dissipation function [N/m 2 s] 
Frequency [s -1 ] 

Ratio of specific heats 



APPENDIX B 


NOZZLE DESIGN 

The nozzle used in this study had a length, b, of 350 milimeters. The domain 
used in the simulations was the top half the a cross-section of the nozzle cut by 
the plane 9 — 0, see figure Bl. The equation representing the nozzle shape was 
actually four separate equations which represent a section of the nozzle. Starting 
from z = 0, the sections are described by a straight fine, a circle, a cubic spline, 
and another straight fine. The location of the divisions are the points al, a2, and 
a3 [figure Bl]. 

The equation of the nozzle is the equation of the radius, r, in meters, which is 
a function of the axial distance z, in milimeters. The equations for each section are 
as follows: 

0 < z < ai, r 
ai < z < a 2 , r 

02 < z < 03, r 

<8 

03 < z < b, r 
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where 


to* = — tan 60° 


m 2 = tan 6 
0 =4° 


ao =0.5 mm 


i2 =2.0 mm 


a 1 


02 


03 


04 


5 — a 0 — R + VR 2 ~3 

Vs 

=ai + V3R/2 

=02 + 2Rm 2 

ARm\ 

=O0 + T" i 


ro =5.0 mm. 
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Figure Bl. Nozzle Design 
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